<!DOCTYPE html>

<html>

<head>

<meta charset="utf-8" />
<meta name="generator" content="pandoc" />
<meta http-equiv="X-UA-Compatible" content="IE=EDGE" />

<meta name="viewport" content="width=device-width, initial-scale=1" />

<meta name="author" content="Edith Ross, Kerstin Haase" />

<meta name="date" content="2020-04-06" />

<title>Multi-sample analysis with ASCAT</title>



<style type="text/css">code{white-space: pre;}</style>
<style type="text/css" data-origin="pandoc">
div.sourceCode { overflow-x: auto; }
table.sourceCode, tr.sourceCode, td.lineNumbers, td.sourceCode {
  margin: 0; padding: 0; vertical-align: baseline; border: none; }
table.sourceCode { width: 100%; line-height: 100%; }
td.lineNumbers { text-align: right; padding-right: 4px; padding-left: 4px; color: #aaaaaa; border-right: 1px solid #aaaaaa; }
td.sourceCode { padding-left: 5px; }
code > span.kw { color: #007020; font-weight: bold; } /* Keyword */
code > span.dt { color: #902000; } /* DataType */
code > span.dv { color: #40a070; } /* DecVal */
code > span.bn { color: #40a070; } /* BaseN */
code > span.fl { color: #40a070; } /* Float */
code > span.ch { color: #4070a0; } /* Char */
code > span.st { color: #4070a0; } /* String */
code > span.co { color: #60a0b0; font-style: italic; } /* Comment */
code > span.ot { color: #007020; } /* Other */
code > span.al { color: #ff0000; font-weight: bold; } /* Alert */
code > span.fu { color: #06287e; } /* Function */
code > span.er { color: #ff0000; font-weight: bold; } /* Error */
code > span.wa { color: #60a0b0; font-weight: bold; font-style: italic; } /* Warning */
code > span.cn { color: #880000; } /* Constant */
code > span.sc { color: #4070a0; } /* SpecialChar */
code > span.vs { color: #4070a0; } /* VerbatimString */
code > span.ss { color: #bb6688; } /* SpecialString */
code > span.im { } /* Import */
code > span.va { color: #19177c; } /* Variable */
code > span.cf { color: #007020; font-weight: bold; } /* ControlFlow */
code > span.op { color: #666666; } /* Operator */
code > span.bu { } /* BuiltIn */
code > span.ex { } /* Extension */
code > span.pp { color: #bc7a00; } /* Preprocessor */
code > span.at { color: #7d9029; } /* Attribute */
code > span.do { color: #ba2121; font-style: italic; } /* Documentation */
code > span.an { color: #60a0b0; font-weight: bold; font-style: italic; } /* Annotation */
code > span.cv { color: #60a0b0; font-weight: bold; font-style: italic; } /* CommentVar */
code > span.in { color: #60a0b0; font-weight: bold; font-style: italic; } /* Information */

</style>
<script>
// apply pandoc div.sourceCode style to pre.sourceCode instead
(function() {
  var sheets = document.styleSheets;
  for (var i = 0; i < sheets.length; i++) {
    if (sheets[i].ownerNode.dataset["origin"] !== "pandoc") continue;
    try { var rules = sheets[i].cssRules; } catch (e) { continue; }
    for (var j = 0; j < rules.length; j++) {
      var rule = rules[j];
      // check if there is a div.sourceCode rule
      if (rule.type !== rule.STYLE_RULE || rule.selectorText !== "div.sourceCode") continue;
      var style = rule.style.cssText;
      // check if color or background-color is set
      if (rule.style.color === '' && rule.style.backgroundColor === '') continue;
      // replace div.sourceCode by a pre.sourceCode rule
      sheets[i].deleteRule(j);
      sheets[i].insertRule('pre.sourceCode{' + style + '}', j);
    }
  }
})();
</script>



<link href="data:text/css;charset=utf-8,body%20%7B%0Abackground%2Dcolor%3A%20%23fff%3B%0Amargin%3A%201em%20auto%3B%0Amax%2Dwidth%3A%20700px%3B%0Aoverflow%3A%20visible%3B%0Apadding%2Dleft%3A%202em%3B%0Apadding%2Dright%3A%202em%3B%0Afont%2Dfamily%3A%20%22Open%20Sans%22%2C%20%22Helvetica%20Neue%22%2C%20Helvetica%2C%20Arial%2C%20sans%2Dserif%3B%0Afont%2Dsize%3A%2014px%3B%0Aline%2Dheight%3A%201%2E35%3B%0A%7D%0A%23TOC%20%7B%0Aclear%3A%20both%3B%0Amargin%3A%200%200%2010px%2010px%3B%0Apadding%3A%204px%3B%0Awidth%3A%20400px%3B%0Aborder%3A%201px%20solid%20%23CCCCCC%3B%0Aborder%2Dradius%3A%205px%3B%0Abackground%2Dcolor%3A%20%23f6f6f6%3B%0Afont%2Dsize%3A%2013px%3B%0Aline%2Dheight%3A%201%2E3%3B%0A%7D%0A%23TOC%20%2Etoctitle%20%7B%0Afont%2Dweight%3A%20bold%3B%0Afont%2Dsize%3A%2015px%3B%0Amargin%2Dleft%3A%205px%3B%0A%7D%0A%23TOC%20ul%20%7B%0Apadding%2Dleft%3A%2040px%3B%0Amargin%2Dleft%3A%20%2D1%2E5em%3B%0Amargin%2Dtop%3A%205px%3B%0Amargin%2Dbottom%3A%205px%3B%0A%7D%0A%23TOC%20ul%20ul%20%7B%0Amargin%2Dleft%3A%20%2D2em%3B%0A%7D%0A%23TOC%20li%20%7B%0Aline%2Dheight%3A%2016px%3B%0A%7D%0Atable%20%7B%0Amargin%3A%201em%20auto%3B%0Aborder%2Dwidth%3A%201px%3B%0Aborder%2Dcolor%3A%20%23DDDDDD%3B%0Aborder%2Dstyle%3A%20outset%3B%0Aborder%2Dcollapse%3A%20collapse%3B%0A%7D%0Atable%20th%20%7B%0Aborder%2Dwidth%3A%202px%3B%0Apadding%3A%205px%3B%0Aborder%2Dstyle%3A%20inset%3B%0A%7D%0Atable%20td%20%7B%0Aborder%2Dwidth%3A%201px%3B%0Aborder%2Dstyle%3A%20inset%3B%0Aline%2Dheight%3A%2018px%3B%0Apadding%3A%205px%205px%3B%0A%7D%0Atable%2C%20table%20th%2C%20table%20td%20%7B%0Aborder%2Dleft%2Dstyle%3A%20none%3B%0Aborder%2Dright%2Dstyle%3A%20none%3B%0A%7D%0Atable%20thead%2C%20table%20tr%2Eeven%20%7B%0Abackground%2Dcolor%3A%20%23f7f7f7%3B%0A%7D%0Ap%20%7B%0Amargin%3A%200%2E5em%200%3B%0A%7D%0Ablockquote%20%7B%0Abackground%2Dcolor%3A%20%23f6f6f6%3B%0Apadding%3A%200%2E25em%200%2E75em%3B%0A%7D%0Ahr%20%7B%0Aborder%2Dstyle%3A%20solid%3B%0Aborder%3A%20none%3B%0Aborder%2Dtop%3A%201px%20solid%20%23777%3B%0Amargin%3A%2028px%200%3B%0A%7D%0Adl%20%7B%0Amargin%2Dleft%3A%200%3B%0A%7D%0Adl%20dd%20%7B%0Amargin%2Dbottom%3A%2013px%3B%0Amargin%2Dleft%3A%2013px%3B%0A%7D%0Adl%20dt%20%7B%0Afont%2Dweight%3A%20bold%3B%0A%7D%0Aul%20%7B%0Amargin%2Dtop%3A%200%3B%0A%7D%0Aul%20li%20%7B%0Alist%2Dstyle%3A%20circle%20outside%3B%0A%7D%0Aul%20ul%20%7B%0Amargin%2Dbottom%3A%200%3B%0A%7D%0Apre%2C%20code%20%7B%0Abackground%2Dcolor%3A%20%23f7f7f7%3B%0Aborder%2Dradius%3A%203px%3B%0Acolor%3A%20%23333%3B%0Awhite%2Dspace%3A%20pre%2Dwrap%3B%20%0A%7D%0Apre%20%7B%0Aborder%2Dradius%3A%203px%3B%0Amargin%3A%205px%200px%2010px%200px%3B%0Apadding%3A%2010px%3B%0A%7D%0Apre%3Anot%28%5Bclass%5D%29%20%7B%0Abackground%2Dcolor%3A%20%23f7f7f7%3B%0A%7D%0Acode%20%7B%0Afont%2Dfamily%3A%20Consolas%2C%20Monaco%2C%20%27Courier%20New%27%2C%20monospace%3B%0Afont%2Dsize%3A%2085%25%3B%0A%7D%0Ap%20%3E%20code%2C%20li%20%3E%20code%20%7B%0Apadding%3A%202px%200px%3B%0A%7D%0Adiv%2Efigure%20%7B%0Atext%2Dalign%3A%20center%3B%0A%7D%0Aimg%20%7B%0Abackground%2Dcolor%3A%20%23FFFFFF%3B%0Apadding%3A%202px%3B%0Aborder%3A%201px%20solid%20%23DDDDDD%3B%0Aborder%2Dradius%3A%203px%3B%0Aborder%3A%201px%20solid%20%23CCCCCC%3B%0Amargin%3A%200%205px%3B%0A%7D%0Ah1%20%7B%0Amargin%2Dtop%3A%200%3B%0Afont%2Dsize%3A%2035px%3B%0Aline%2Dheight%3A%2040px%3B%0A%7D%0Ah2%20%7B%0Aborder%2Dbottom%3A%204px%20solid%20%23f7f7f7%3B%0Apadding%2Dtop%3A%2010px%3B%0Apadding%2Dbottom%3A%202px%3B%0Afont%2Dsize%3A%20145%25%3B%0A%7D%0Ah3%20%7B%0Aborder%2Dbottom%3A%202px%20solid%20%23f7f7f7%3B%0Apadding%2Dtop%3A%2010px%3B%0Afont%2Dsize%3A%20120%25%3B%0A%7D%0Ah4%20%7B%0Aborder%2Dbottom%3A%201px%20solid%20%23f7f7f7%3B%0Amargin%2Dleft%3A%208px%3B%0Afont%2Dsize%3A%20105%25%3B%0A%7D%0Ah5%2C%20h6%20%7B%0Aborder%2Dbottom%3A%201px%20solid%20%23ccc%3B%0Afont%2Dsize%3A%20105%25%3B%0A%7D%0Aa%20%7B%0Acolor%3A%20%230033dd%3B%0Atext%2Ddecoration%3A%20none%3B%0A%7D%0Aa%3Ahover%20%7B%0Acolor%3A%20%236666ff%3B%20%7D%0Aa%3Avisited%20%7B%0Acolor%3A%20%23800080%3B%20%7D%0Aa%3Avisited%3Ahover%20%7B%0Acolor%3A%20%23BB00BB%3B%20%7D%0Aa%5Bhref%5E%3D%22http%3A%22%5D%20%7B%0Atext%2Ddecoration%3A%20underline%3B%20%7D%0Aa%5Bhref%5E%3D%22https%3A%22%5D%20%7B%0Atext%2Ddecoration%3A%20underline%3B%20%7D%0A%0Acode%20%3E%20span%2Ekw%20%7B%20color%3A%20%23555%3B%20font%2Dweight%3A%20bold%3B%20%7D%20%0Acode%20%3E%20span%2Edt%20%7B%20color%3A%20%23902000%3B%20%7D%20%0Acode%20%3E%20span%2Edv%20%7B%20color%3A%20%2340a070%3B%20%7D%20%0Acode%20%3E%20span%2Ebn%20%7B%20color%3A%20%23d14%3B%20%7D%20%0Acode%20%3E%20span%2Efl%20%7B%20color%3A%20%23d14%3B%20%7D%20%0Acode%20%3E%20span%2Ech%20%7B%20color%3A%20%23d14%3B%20%7D%20%0Acode%20%3E%20span%2Est%20%7B%20color%3A%20%23d14%3B%20%7D%20%0Acode%20%3E%20span%2Eco%20%7B%20color%3A%20%23888888%3B%20font%2Dstyle%3A%20italic%3B%20%7D%20%0Acode%20%3E%20span%2Eot%20%7B%20color%3A%20%23007020%3B%20%7D%20%0Acode%20%3E%20span%2Eal%20%7B%20color%3A%20%23ff0000%3B%20font%2Dweight%3A%20bold%3B%20%7D%20%0Acode%20%3E%20span%2Efu%20%7B%20color%3A%20%23900%3B%20font%2Dweight%3A%20bold%3B%20%7D%20%0Acode%20%3E%20span%2Eer%20%7B%20color%3A%20%23a61717%3B%20background%2Dcolor%3A%20%23e3d2d2%3B%20%7D%20%0A" rel="stylesheet" type="text/css" />




</head>

<body>




<h1 class="title toc-ignore">Multi-sample analysis with ASCAT</h1>
<h4 class="author">Edith Ross, Kerstin Haase</h4>
<h4 class="date">2020-04-06</h4>



<div id="vignette-info" class="section level2">
<h2>Vignette Info</h2>
<p>This vignette demonstrates how to use ASCAT to analyse multiple phylogenetically related samples. For the general usage of ASCAT including parameters that are not specific to multi-sample analysis please refer to the <a href="https://www.crick.ac.uk/peter-van-loo/software/ASCAT">ASCAT webpage</a> and the <a href="https://github.com/Crick-CancerGenomics/ascat/blob/master/ExampleData/ASCAT_examplePipeline.R">example pipeline</a>.</p>
</div>
<div id="using-ascat-for-multi-sample-analyses" class="section level2">
<h2>Using ASCAT for multi-sample analyses</h2>
<p>We start by loading the ASCAT package.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(ASCAT)</code></pre></div>
<p>Next we load the data.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">ascat.bcMulti &lt;-<span class="st"> </span><span class="kw">ascat.loadData</span>(
  <span class="dt">Tumor_LogR_file =</span> <span class="kw">system.file</span>(<span class="st">&quot;extdata&quot;</span>, <span class="st">&quot;tumour.logR.txt&quot;</span>, <span class="dt">package=</span><span class="st">&quot;ASCAT&quot;</span>),
  <span class="dt">Tumor_BAF_file =</span> <span class="kw">system.file</span>(<span class="st">&quot;extdata&quot;</span>, <span class="st">&quot;tumour.BAF.txt&quot;</span>, <span class="dt">package=</span><span class="st">&quot;ASCAT&quot;</span>),
  <span class="dt">Germline_LogR_file =</span> <span class="kw">system.file</span>(<span class="st">&quot;extdata&quot;</span>, <span class="st">&quot;singlenormal.logR.txt&quot;</span>, <span class="dt">package=</span><span class="st">&quot;ASCAT&quot;</span>),
  <span class="dt">Germline_BAF_file =</span> <span class="kw">system.file</span>(<span class="st">&quot;extdata&quot;</span>, <span class="st">&quot;singlenormal.BAF.txt&quot;</span>, <span class="dt">package=</span><span class="st">&quot;ASCAT&quot;</span>))</code></pre></div>
<pre><code>## [1] Reading Tumor LogR data...
## [1] Reading Tumor BAF data...
## [1] Reading Germline LogR data...
## [1] Reading Germline BAF data...
## [1] Registering SNP locations...
## [1] Splitting genome in distinct chunks...</code></pre>
<p>Both <tt>Tumor_LogR_file</tt> and <tt>Tumor_BAF_file</tt> are expected to contain a column for each of the samples to analyse.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">head</span>(ascat.bcMulti<span class="op">$</span>Tumor_LogR)                          </code></pre></div>
<pre><code>##            S1       S2
## SNP1  0.03615 -1.03950
## SNP2  0.14998 -0.79433
## SNP3 -0.00891 -0.76137
## SNP4  0.40188 -0.67521
## SNP5  0.14902 -0.72980
## SNP6  0.24118 -1.11302</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">head</span>(ascat.bcMulti<span class="op">$</span>Tumor_BAF)       </code></pre></div>
<pre><code>##           S1      S2
## SNP1 0.51596 0.99262
## SNP2 0.67903 0.00255
## SNP3 1.00000 1.00000
## SNP4 0.00000 0.00000
## SNP5 1.00000 1.00000
## SNP6 0.45572 0.04925</code></pre>
<p>The next step is to run the segmentation. When analysing phylogenetically related samples, it is expected that some of the copy number segment boundaries are shared between samples. In this case a joint segmentation of all samples is recommended. The synthetic data set used in this example was also simulated with partly shared segment boundaries. The ground truth copy number plots of the two samples we are going to analyse are shown in the following plots.</p>
<p><img src="" style="display: block; margin: auto;" /><img src="" style="display: block; margin: auto;" /></p>
<p>The multi-sample segmentation algorithm can be run using the function <tt>ascat.asmultipcf</tt>.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">ascat.bcMulti &lt;-<span class="st"> </span><span class="kw">ascat.asmultipcf</span>(ascat.bcMulti,<span class="dt">penalty =</span> <span class="dv">5</span>)</code></pre></div>
<pre><code>## [1] &quot;Segmentlength 5&quot;</code></pre>
<p>Finally ASCAT can be run on the segmented data set.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">ascat.outputMulti =<span class="st"> </span><span class="kw">ascat.runAscat</span>(ascat.bcMulti)</code></pre></div>
</div>
<div id="comparison-with-single-sample-segmentation" class="section level2">
<h2>Comparison with single sample segmentation</h2>
<p>Finally, we compare our result to that of standard single sample segmentation using <tt>ascat.aspcf</tt>.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">ascat.bc =<span class="st"> </span><span class="kw">ascat.loadData</span>(<span class="kw">system.file</span>(<span class="st">&quot;extdata&quot;</span>, <span class="st">&quot;tumour.logR.txt&quot;</span>, <span class="dt">package=</span><span class="st">&quot;ASCAT&quot;</span>),
                          <span class="kw">system.file</span>(<span class="st">&quot;extdata&quot;</span>, <span class="st">&quot;tumour.BAF.txt&quot;</span>, <span class="dt">package=</span><span class="st">&quot;ASCAT&quot;</span>),
                          <span class="kw">system.file</span>(<span class="st">&quot;extdata&quot;</span>, <span class="st">&quot;normal.logR.txt&quot;</span>, <span class="dt">package=</span><span class="st">&quot;ASCAT&quot;</span>),
                          <span class="kw">system.file</span>(<span class="st">&quot;extdata&quot;</span>, <span class="st">&quot;normal.BAF.txt&quot;</span>, <span class="dt">package=</span><span class="st">&quot;ASCAT&quot;</span>))</code></pre></div>
<pre><code>## [1] Reading Tumor LogR data...
## [1] Reading Tumor BAF data...
## [1] Reading Germline LogR data...
## [1] Reading Germline BAF data...
## [1] Registering SNP locations...
## [1] Splitting genome in distinct chunks...</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">ascat.bc =<span class="st"> </span><span class="kw">ascat.aspcf</span>(ascat.bc,<span class="dt">penalty =</span> <span class="dv">25</span>)</code></pre></div>
<pre><code>## [1] Sample S1 (1/2)
## [1] Sample S2 (2/2)</code></pre>
<p>Note that in the single-sample case the same segmentation sensitivity is achieved with a higher penalty parameter compared to the multi-sample case. This means, when switching from single- to multi-sample segmentation, the penalty parameter needs to be lowered to maintain a similar sensitivity.</p>
<p>We plot the segment boundaries inferred for each of the two samples by multi- and single-sample segmentation.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">plot.segments</span>(<span class="dt">v1=</span><span class="kw">cumsum</span>(<span class="kw">rle</span>(ascat.bc<span class="op">$</span>Tumor_LogR_segmented[,<span class="dv">1</span>])<span class="op">$</span>lengths),
              <span class="dt">v2=</span><span class="kw">cumsum</span>(<span class="kw">rle</span>(ascat.bc<span class="op">$</span>Tumor_LogR_segmented[,<span class="dv">2</span>])<span class="op">$</span>lengths),
              <span class="dt">main=</span><span class="st">&quot;Single-sample segmentation&quot;</span>)

<span class="kw">plot.segments</span>(<span class="dt">v1=</span><span class="kw">cumsum</span>(<span class="kw">rle</span>(ascat.bcMulti<span class="op">$</span>Tumor_LogR_segmented[,<span class="dv">1</span>])<span class="op">$</span>lengths),
              <span class="dt">v2=</span><span class="kw">cumsum</span>(<span class="kw">rle</span>(ascat.bcMulti<span class="op">$</span>Tumor_LogR_segmented[,<span class="dv">2</span>])<span class="op">$</span>lengths),
              <span class="dt">main=</span><span class="st">&quot;Multi-sample segmentation&quot;</span>)</code></pre></div>
<p><img src="" style="display: block; margin: auto;" /><img src="" style="display: block; margin: auto;" /> In case of single-sample segmentation the inferred positions of most of the shared segment boundaries vary slightly between the two samples, whereas the multi-sample segmentation infers a common breakpoint when there is no significant evidence that the boundaries differ between samples.</p>
</div>
<div id="comparison-with-another-multi-sample-copy-number-method" class="section level2">
<h2>Comparison with another multi-sample copy number method</h2>
<p>In order to compare asmultipcf segmentation to other copy number inference methods, we ran ASCAT on two samples from two patients from a metastatic prostate cancer study (<a href="https://www.ncbi.nlm.nih.gov/pubmed/25830880">Gundem et al</a>). For the same study, copy number profiles are available from HATCHet (<a href="https://github.com/raphael-group/hatchet">github repository</a>; <a href="https://www.biorxiv.org/content/10.1101/496174v1">publication</a>).</p>
<p>Loading all neccessary packages.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(ggplot2)
<span class="kw">library</span>(plyr)</code></pre></div>
<p>Because HATCHet can model a mixture of copy number states but ASCAT only detects the major clone, we need to define the major clone in the HATCHet data.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">getMajorClone&lt;-<span class="cf">function</span>(row){
  maj&lt;-<span class="kw">names</span>(<span class="kw">which.max</span>(row[<span class="kw">c</span>(<span class="dv">8</span>,<span class="dv">10</span>,<span class="dv">12</span>)]))
  maj.cn&lt;-row[(<span class="kw">match</span>(maj,<span class="kw">names</span>(row))<span class="op">-</span><span class="dv">1</span>)]
  <span class="kw">return</span>(maj.cn)
}</code></pre></div>
<p>Next we define a function to load HATCHet data, define the major clone for each segment, i.e. row, and split copy number annotation from “a|b” format into two columns.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">readHatchetFile&lt;-<span class="cf">function</span>(path){
  hatchet.raw&lt;-<span class="kw">read.table</span>(path, <span class="dt">header=</span>F, <span class="dt">sep=</span><span class="st">&quot;</span><span class="ch">\t</span><span class="st">&quot;</span>, <span class="dt">stringsAsFactors =</span> F)
  
  ##find major clone
  major.cn&lt;-<span class="kw">apply</span>(hatchet.raw, <span class="dv">1</span>, <span class="cf">function</span>(z) <span class="kw">getMajorClone</span>(z))
  hatchet.raw<span class="op">$</span>major.cn&lt;-major.cn
  
  hatchet.short&lt;-hatchet.raw[<span class="kw">c</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">4</span>,<span class="dv">13</span>)]
  <span class="kw">names</span>(hatchet.short)&lt;-<span class="kw">c</span>(<span class="st">&quot;chr&quot;</span>,<span class="st">&quot;start&quot;</span>,<span class="st">&quot;end&quot;</span>,<span class="st">&quot;sample&quot;</span>,<span class="st">&quot;hatchet.cn&quot;</span>)
  
  major.cn&lt;-<span class="kw">unlist</span>(<span class="kw">lapply</span>(<span class="kw">strsplit</span>(hatchet.short<span class="op">$</span>hatchet.cn, <span class="st">&quot;</span><span class="ch">\\</span><span class="st">|&quot;</span>), <span class="cf">function</span>(z) z[<span class="dv">1</span>]))
  minor.cn&lt;-<span class="kw">unlist</span>(<span class="kw">lapply</span>(<span class="kw">strsplit</span>(hatchet.short<span class="op">$</span>hatchet.cn, <span class="st">&quot;</span><span class="ch">\\</span><span class="st">|&quot;</span>), <span class="cf">function</span>(z) z[<span class="dv">2</span>]))
  
  hatchet.short<span class="op">$</span>major.cn&lt;-major.cn
  hatchet.short<span class="op">$</span>minor.cn&lt;-minor.cn
  
  <span class="kw">return</span>(hatchet.short[<span class="op">-</span><span class="dv">5</span>])
}</code></pre></div>
<p>In order to provide a quantitative measure of how much ASCAT using asmultipcf and HATCHet agree with regard to the sample segmentation, we analysed what fraction of asmultipcf breakpoints has a corresponding HATCHet breakpoint. First, we centred the ASCAT breakpoints in the middle of two segments.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">returnBreakpoints&lt;-<span class="cf">function</span>(df.chr){
  bp.df&lt;-<span class="ot">NULL</span>
  df.sorted&lt;-df.chr[<span class="kw">order</span>(df.chr<span class="op">$</span>start),]
  <span class="cf">for</span>(i <span class="cf">in</span> <span class="dv">1</span><span class="op">:</span>(<span class="kw">nrow</span>(df.sorted)<span class="op">-</span><span class="dv">1</span>)){
    first&lt;-df.sorted[i,]
    second&lt;-df.sorted[(i<span class="op">+</span><span class="dv">1</span>),]
    ##verify real cn bp
    <span class="cf">if</span>(first<span class="op">$</span>major.cn<span class="op">==</span>second<span class="op">$</span>major.cn <span class="op">&amp;</span><span class="st"> </span>first<span class="op">$</span>minor.cn<span class="op">==</span>second<span class="op">$</span>minor.cn){
      <span class="kw">print</span>(<span class="st">&quot;Consecutive rows have identical CN!&quot;</span>)
    }
    <span class="cf">else</span>{
      bp.df&lt;-<span class="kw">rbind</span>(bp.df, <span class="kw">data.frame</span>(<span class="dt">bp.location=</span>(first<span class="op">$</span>end<span class="op">+</span>((second<span class="op">$</span>start<span class="op">-</span>first<span class="op">$</span>end)<span class="op">/</span><span class="dv">2</span>))))
    }
  }
  <span class="kw">return</span>(bp.df)
}</code></pre></div>
<p>HATCHet output is provided for 50kb bins and not specifically segmented. Hence, in order to compare consecutive segments with the same major clonal copy number, we are merging neighbouring segments with matching allele specific copy number values. Then we calculate the distance of the closest HATCHet breakpoint for every asmultipcf breakpoint.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">findMatchingBPs&lt;-<span class="cf">function</span>(ascat.bps, hatchet.segs){
  ##limit hatchet df to true CN segs
  hatchet.true&lt;-<span class="ot">NULL</span>
  hatchet.sorted&lt;-hatchet.segs[<span class="kw">order</span>(hatchet.segs<span class="op">$</span>start),]
  
  rles&lt;-<span class="kw">rle</span>(<span class="kw">paste0</span>(hatchet.sorted<span class="op">$</span>major.cn,<span class="st">&quot;_&quot;</span>,hatchet.sorted<span class="op">$</span>minor.cn))<span class="op">$</span>lengths
  idx&lt;-<span class="dv">1</span>
  <span class="cf">for</span>(k <span class="cf">in</span> rles){
    hatchet.true&lt;-<span class="kw">rbind</span>(hatchet.true, <span class="kw">data.frame</span>(<span class="dt">start=</span>hatchet.sorted[idx, <span class="st">&quot;start&quot;</span>], <span class="dt">end=</span>hatchet.sorted[(idx<span class="op">+</span>k<span class="op">-</span><span class="dv">1</span>), <span class="st">&quot;end&quot;</span>]))
    idx&lt;-(idx<span class="op">+</span>k)
  }
  
  ascat.flag&lt;-<span class="ot">NULL</span>
  <span class="cf">for</span>(i <span class="cf">in</span> <span class="dv">1</span><span class="op">:</span><span class="kw">nrow</span>(ascat.bps)){
    ascat.bp&lt;-ascat.bps[i,<span class="st">&quot;bp.location&quot;</span>]
    min.dist&lt;-<span class="kw">min</span>(<span class="kw">c</span>(<span class="kw">abs</span>(ascat.bp<span class="op">-</span>hatchet.true<span class="op">$</span>end),<span class="kw">abs</span>(ascat.bp<span class="op">-</span>hatchet.true<span class="op">$</span>start)))
    
    ascat.flag&lt;-<span class="kw">rbind</span>(ascat.flag, <span class="kw">data.frame</span>(<span class="dt">bkp=</span>ascat.bp, <span class="dt">min.dist=</span>min.dist))
  }
  <span class="kw">return</span>(ascat.flag)
}</code></pre></div>
<p>This is the main function visualising the copy number of asmultipcf and HATCHet for all four analysed samples and calculates the fraction of ASCAT breakpoints with a HATCHet breakpoint fewer than 50kb bases away (size of HATCHet bins).</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">compareCN&lt;-<span class="cf">function</span>(cn.ascat, cn.hatchet){
  ##split into samples
  shared.samples&lt;-<span class="kw">intersect</span>(cn.ascat<span class="op">$</span>sample, cn.hatchet<span class="op">$</span>sample)
  
  <span class="cf">for</span>(s <span class="cf">in</span> shared.samples){
    ascat.sub&lt;-cn.ascat[cn.ascat<span class="op">$</span>sample<span class="op">==</span>s,]
    hatchet.sub&lt;-cn.hatchet[cn.hatchet<span class="op">$</span>sample<span class="op">==</span>s,]
    ##move values slightly for plotting
    ascat.plot&lt;-<span class="st"> </span>ascat.sub
    hatchet.plot&lt;-<span class="st"> </span>hatchet.sub
    ascat.plot<span class="op">$</span>major.cn&lt;-(<span class="kw">as.numeric</span>(ascat.plot<span class="op">$</span>major.cn)<span class="op">+</span><span class="fl">0.1</span>)
    hatchet.plot<span class="op">$</span>major.cn&lt;-(<span class="kw">as.numeric</span>(hatchet.plot<span class="op">$</span>major.cn)<span class="op">-</span><span class="fl">0.1</span>)
    
    ##visualise
    joint.df&lt;-<span class="kw">rbind</span>(ascat.plot,hatchet.plot)
    joint.df<span class="op">$</span>chr&lt;-<span class="kw">factor</span>(joint.df<span class="op">$</span>chr, <span class="dt">levels=</span><span class="kw">c</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">22</span>,<span class="st">&quot;X&quot;</span>,<span class="st">&quot;Y&quot;</span>))
    p1&lt;-<span class="kw">ggplot</span>(joint.df, <span class="kw">aes</span>(<span class="dt">x=</span>start, <span class="dt">y=</span>(<span class="kw">as.numeric</span>(major.cn)<span class="op">+</span><span class="kw">as.numeric</span>(minor.cn)), <span class="dt">xend=</span>end, <span class="dt">yend=</span>(<span class="kw">as.numeric</span>(major.cn)<span class="op">+</span><span class="kw">as.numeric</span>(minor.cn)), <span class="dt">col=</span><span class="kw">as.factor</span>(method)))<span class="op">+</span>
<span class="st">      </span><span class="kw">geom_segment</span>()<span class="op">+</span>
<span class="st">      </span><span class="kw">theme</span>(<span class="dt">axis.text.x=</span><span class="kw">element_text</span>(<span class="dt">angle=</span><span class="dv">45</span>,<span class="dt">hjust=</span><span class="dv">1</span>), <span class="dt">panel.background =</span> <span class="kw">element_blank</span>(), <span class="dt">panel.grid.major =</span> <span class="kw">element_line</span>(<span class="dt">color=</span><span class="st">&quot;grey80&quot;</span>))<span class="op">+</span>
<span class="st">      </span><span class="kw">ggtitle</span>(<span class="kw">paste0</span>(<span class="st">&quot;Sample &quot;</span>,s))<span class="op">+</span>
<span class="st">      </span><span class="kw">scale_color_manual</span>(<span class="dt">values=</span><span class="kw">c</span>(<span class="st">&quot;cyan4&quot;</span>,<span class="st">&quot;coral3&quot;</span>), <span class="dt">name=</span><span class="st">&quot;&quot;</span>)<span class="op">+</span>
<span class="st">      </span><span class="kw">scale_x_continuous</span>(<span class="dt">name=</span><span class="st">&quot;&quot;</span>)<span class="op">+</span>
<span class="st">      </span><span class="kw">scale_y_continuous</span>(<span class="dt">name=</span><span class="st">&quot;Copy number&quot;</span>, <span class="dt">limits=</span>(<span class="kw">c</span>(<span class="dv">0</span>,<span class="dv">15</span>)))<span class="op">+</span>
<span class="st">      </span><span class="kw">facet_wrap</span>(.<span class="op">~</span>chr, <span class="dt">scales =</span> <span class="st">&quot;free_x&quot;</span>, <span class="dt">nrow =</span> <span class="dv">3</span>)
    <span class="kw">print</span>(p1)
    
    
    ##give quality in measure: what fraction of ASCAT bp has hatchet bp in &lt; threshold vincinity?
    bps&lt;-<span class="kw">ddply</span>(ascat.sub, .(chr), <span class="cf">function</span>(z) <span class="kw">returnBreakpoints</span>(z))
    bp.list&lt;-<span class="kw">split</span>(bps, bps<span class="op">$</span>chr)
    
    bkp.list&lt;-<span class="ot">NULL</span>
    <span class="cf">for</span>(c <span class="cf">in</span> <span class="kw">names</span>(bp.list)){
      <span class="cf">if</span>(c<span class="op">!=</span><span class="st">&quot;X&quot;</span>){
        bkp.list&lt;-<span class="kw">rbind</span>(bkp.list,<span class="kw">findMatchingBPs</span>(bp.list[[c]], hatchet.sub[hatchet.sub<span class="op">$</span>chr<span class="op">==</span>c,]))
      }
    }
    ##use bin size of HATCHet as threshold
    threshold=<span class="dv">50000</span>
    bkp.list<span class="op">$</span>match&lt;-<span class="kw">ifelse</span>(bkp.list<span class="op">$</span>min.dist<span class="op">&lt;</span>threshold,T,F)
    <span class="kw">print</span>(s)
    <span class="kw">print</span>(<span class="kw">prop.table</span>(<span class="kw">table</span>(bkp.list<span class="op">$</span>match)))
  }
}</code></pre></div>
<p>Now we can carry out the comparison. Load the HATCHet segmentation for patient A32.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">path&lt;-<span class="st">&quot;https://raw.githubusercontent.com/raphael-group/hatchet-paper/master/cancer/prostate/A32/best.seg.ucn&quot;</span>
hatchet.seg&lt;-<span class="kw">readHatchetFile</span>(path)
hatchet.seg<span class="op">$</span>method&lt;-<span class="st">&quot;HATCHet&quot;</span>
hatchet.seg<span class="op">$</span>sample&lt;-<span class="kw">gsub</span>(<span class="st">&quot;-&quot;</span>,<span class="st">&quot;&quot;</span>,hatchet.seg<span class="op">$</span>sample)</code></pre></div>
<p>Now we can load the asmultipcf results for patient A32.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">ascat.seg&lt;-<span class="kw">read.table</span>(<span class="kw">system.file</span>(<span class="st">&quot;extdata&quot;</span>, <span class="st">&quot;A32.fast.segments.txt&quot;</span>, <span class="dt">package=</span><span class="st">&quot;ASCAT&quot;</span>), <span class="dt">header=</span>T, <span class="dt">sep=</span><span class="st">&quot;</span><span class="ch">\t</span><span class="st">&quot;</span>, <span class="dt">stringsAsFactors =</span> F)
ascat.seg<span class="op">$</span>tmp&lt;-<span class="kw">unlist</span>(<span class="kw">lapply</span>(<span class="kw">strsplit</span>(ascat.seg<span class="op">$</span>sample, <span class="st">&quot;-&quot;</span>), <span class="cf">function</span>(z) z[<span class="dv">1</span>]))
ascat.seg&lt;-ascat.seg[<span class="op">-</span><span class="dv">1</span>]
<span class="kw">names</span>(ascat.seg)&lt;-<span class="kw">c</span>(<span class="st">&quot;chr&quot;</span>,<span class="st">&quot;start&quot;</span>,<span class="st">&quot;end&quot;</span>,<span class="st">&quot;major.cn&quot;</span>,<span class="st">&quot;minor.cn&quot;</span>,<span class="st">&quot;sample&quot;</span>)
ascat.seg<span class="op">$</span>method&lt;-<span class="st">&quot;ASCAT&quot;</span></code></pre></div>
<p>Finally, we can create the visual and segmentation comparison between the copy number of asmultipcf and HATCHet.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">compareCN</span>(ascat.seg, hatchet.seg)</code></pre></div>
<p><img src="" style="display: block; margin: auto;" /></p>
<pre><code>## [1] &quot;A32C&quot;
## 
##     FALSE      TRUE 
## 0.3591331 0.6408669</code></pre>
<p><img src="" style="display: block; margin: auto;" /></p>
<pre><code>## [1] &quot;A32A&quot;
## 
##     FALSE      TRUE 
## 0.4606061 0.5393939</code></pre>
<p>So for both samples of patient A32, between 54% and 64% of asmultipcf breakpoints have a corresponding HATCHet breakpoint.</p>
<p>We can now repeat the same analysis for patient A17.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">path&lt;-<span class="st">&quot;https://raw.githubusercontent.com/raphael-group/hatchet-paper/master/cancer/prostate/A17/best.seg.ucn&quot;</span>
hatchet.seg&lt;-<span class="kw">readHatchetFile</span>(path)
hatchet.seg<span class="op">$</span>method&lt;-<span class="st">&quot;HATCHet&quot;</span>

ascat.seg&lt;-<span class="kw">read.table</span>(<span class="kw">system.file</span>(<span class="st">&quot;extdata&quot;</span>, <span class="st">&quot;A17.fast.segments.txt&quot;</span>, <span class="dt">package=</span><span class="st">&quot;ASCAT&quot;</span>), <span class="dt">header=</span>T, <span class="dt">sep=</span><span class="st">&quot;</span><span class="ch">\t</span><span class="st">&quot;</span>, <span class="dt">stringsAsFactors =</span> F)
ascat.seg<span class="op">$</span>tmp&lt;-<span class="kw">unlist</span>(<span class="kw">lapply</span>(<span class="kw">strsplit</span>(ascat.seg<span class="op">$</span>sample, <span class="st">&quot;-&quot;</span>), <span class="cf">function</span>(z) z[<span class="dv">1</span>]))
ascat.seg&lt;-ascat.seg[<span class="op">-</span><span class="dv">1</span>]
<span class="kw">names</span>(ascat.seg)&lt;-<span class="kw">c</span>(<span class="st">&quot;chr&quot;</span>,<span class="st">&quot;start&quot;</span>,<span class="st">&quot;end&quot;</span>,<span class="st">&quot;major.cn&quot;</span>,<span class="st">&quot;minor.cn&quot;</span>,<span class="st">&quot;sample&quot;</span>)
ascat.seg<span class="op">$</span>method&lt;-<span class="st">&quot;ASCAT&quot;</span></code></pre></div>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">compareCN</span>(ascat.seg, hatchet.seg)</code></pre></div>
<p><img src="" style="display: block; margin: auto;" /></p>
<pre><code>## [1] &quot;A17F&quot;
## 
##     FALSE      TRUE 
## 0.2039801 0.7960199</code></pre>
<p><img src="" style="display: block; margin: auto;" /></p>
<pre><code>## [1] &quot;A17A&quot;
## 
##     FALSE      TRUE 
## 0.1955017 0.8044983</code></pre>
<p>For both samples of patient A17, around 80% of asmultipcf breakpoints have a corresponding HATCHet breakpoint.</p>
</div>



<!-- code folding -->


<!-- dynamically load mathjax for compatibility with self-contained -->
<script>
  (function () {
    var script = document.createElement("script");
    script.type = "text/javascript";
    script.src  = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML";
    document.getElementsByTagName("head")[0].appendChild(script);
  })();
</script>

</body>
</html>
